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Abstract 

The 5 = 1 antiferromagnetic Heisenberg chain both with and without single 
ion anisotropy is studied. Using the recently proposed density matrix renor- 
malization group technique we calculate the energy gaps as well as several 
different correlation functions. The two gaps, A||, A_|_, along with associated 
correlation lengths and velocities are determined. The numerical results are 
shown to be in good agreement with theoretical predictions derived from the 
nonlinear sigma model and a free boson model. We also study the 5 = 1/2 
excitations that occur at the ends of open chains; in particular we study the 
behavior associated with open boundary conditions, using a model of 5 = 1/2 
spins coupled to the free bosons. 
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I. INTRODUCTION 



Following the initial conjecture of HaldaneB it is by now well established both experimen- 
tallyi'i and theoreticallyoi that integer spin antiferromagnetic chains have a gap as opposed 
to chains with half-integer spin. The spectrum is that of a disordered singlet ground-state 
with a gap, A, to a triplet excitation. The disordered Haldane phase is known to have a 
hidden order parameter, the so called string order parameter^. One of the best candidates 
for such a Haldane gap system is Ni(C2H8N2)2N02(C104) (NENP). The ratio of inter- to 
intra chain coupling in this compound is estimatedO to be less than 4 x 10"'* and no ordering 
has so far been observed at any accessible temperatures. NENP is not a completely isotropic 
compound. As revealed by neutron scattering experiments the magnon spectrum around 
wave vector vr is split into three branchesilll due to the effects of anisotropics. Electron 
spm resonance studieMl have observed transitions between the different components of 
the triplet confirming the expected nature of the spectrum. Thus, the isotropic Haldane 
gap. A, is effectively split in three. Single ion anisotropy, D{S-)^, splits the triplet into two 
modes for fluctuations perpendicular and parallel to the chain, respectively. In addition, 
in-plane anisotropy, of the form E{{SfY ~ (Sf)'^), splits the perpendicular mode by a small 
amount. At present the best estimates of the gaps in NENPll3ilil are 2.50, 1.25, and 1.05 
meV, respectively. The dispersion relation and dynamical structure factor have also been 
measured for NENF^'EHll. The largest of the anisotropics, and the only one we shall be 
concerned with here, is the single-ion anisotropy, D. A simple Hamiltonian capturing the 
essential physics including the single ion anisotropy is then 

H = JJ2{SrS^+l + D{S!f}. (1.1) 

i 

As a function of D the lowest lying gap will decrease. At a critical value Dc ^ 1 the gap 
closes and a phase transition between the disordered Haldane phase and a "large-D" phase 
occur^i. For a discussion of the phase diagram see Ref. |10[ By fitting exact diagonalization 
results on chains of length up to 18 to inelastic neutron scattering results Golinelli et allll 
have determined the optimal values J = 3.75 meV and D/ J = 0.18. This value of D is then 
clearly in the disordered Haldane phase. Even including inter chain coupling, which changes 
the value of Dc, D / J = 0.18 remains in the disordered phased, consistent with the fact that 
no ordering is observed in NENP at finite temperatures. Recently it has been shown that 
the Hamiltonian Eq. (|TTT|) should be corrected to account for the staggered structure of the 
NENP chainsS, a fact that also leads to a staggered off-diagonal part in the gyromagnetic 
tensoiii. 

The longest chains one has studied by exact diagonalization have 18 siteSiiiS. While 
it is possible to study somewhat longer chains using quantum Monte Carlo (QMC) tech- 
niquesi'i the statistical errors on these results will limit the detailed interpretation of the 
data. Recently White@ has proposed a density matrix renormalization group (DMRG) 
method that allows one to study significantly longer chains with much higher accuracy than 
what is obtainable using QMC methods. At present it is not known how to work with the 
momentum as a quantum number within the DMRG framework although this would be very 
desirable and one is therefore restricted to study quantities defined in real space. Another 
less important restriction imposed by the DMRG is the fact that it works notably better for 
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chains with open ends. Recent progresaHj indicates that open ended systems can be tuned 
to mimic the behavior of systems subject to periodic boundary conditions. In this paper 
we calculate equal time correlation functions in real space for the S = 1 antiferromagnetic 
chain for systems of up to 100 spins using the DMRG method. The numerical results are 
then compared to predictions obtained from the nonlinear sigma (NLa) model and a free 
boson model. Some of our results have previously been presented elsewhereil. 

Equal time correlation functions have been calculated by QMC techniques&B and by 
exact diagonalizationH for chains of length up to L=18. The lowest energy excitation 
as a function of momentum has been calculated by QMC techniques&i and by exact diago- 
nalizatio nliilil for chains of length up to L=18. The dynamic structure factor has been 
studied by QMC techniqueJUSilil. The dynamical structure factor, including anisotropy, 
for both in plane and out of plane modes has been studied by Golinelli et. ala using exact 
diagonalization on chains with L < 16, and by Haas et. al.@ for the out of plane mode 
with L < 18. While generally applicable the current state of the art exact diagonalization 
studies are limited by the fact that only very short chains can be studied, thus restricting the 
accessible momentum values considerably. Secondly, finite size effects are likely to be large. 
The finite size effects can be improved upon by the application of extrapolation techniques, 
(Shanks transformations), if enough data points are available. Thus, values for the gaps can 
be extracted to a very good precision. However, for the structure factor, this technique can 
only be applied for a few values of the momentum. The QMC results, while applicable to 
considerably longer chains, have statistical errors. The extraction of dynamical properties 
from QMC results requires the analytical continuation of results obtained at imaginary fre- 
quencies to real frequencies. This is typically done using maximum entropy methods and 
is not a trivial matter. While dynamical properties are not accessible yet using the DMRG 
method, it is possible to obtain equal time correlations for very long chains with an unprece- 
dented accuracy. This huge gain in detail and precision allows a detailed comparison to field 
theory and experiment not previously possible. 

We shall here only be concerned with chains of even length subject to open boundary 
conditions. All of our calculations are for chains described by the Hamiltonian Eq. ( p. . 1|) 
with J/kf, =1 and D/J = 0.18 or D/J = 0.0. For these chains the ground-state is a singlet 
with even parity, 0"^. Above the ground-state is an exponentially low-lying triplet, 1~. In 
the thermodynamic limit the triplet and the singlet become degenerate and the ground-state 
four-fold degenerate. The bulk of the calculations of the correlation functions presented here 
will be for a 100 site chain in the 1~ state since this system conceptually is the simplest. 
We have repeated the calculations for the 0"*" state and also performed the same calculations 
for a 60 site system. In none of these cases were the correlation functions seen to differ 
markedly from what we will describe below, aside from boundary effects. In particular, all 
the correlation lengths obtained were consistent. 

We implement the DMRG using density matrices of the size 243 x 243 keeping 81 eigen- 
vectors of these matrices at each iteration. First density matrices representing systems of 
different sizes are computed using an infinite lattice method, then these are combined to form 
a system of a fixed size. L=100, for which all the correlation functions are calculated using 
a finite lattice methodcS. For a discussion of the numerical procedure we refer the reader to 



Ref. 25,p8[ The DMRG method for open chains leaves two good quantum numbers the total 



component, S^, and the parity, P, corresponding to a reflection about the midpoint of 
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the chain. These are conserved under iteration and it is therefore possible to work within 
a subspace defined by these two quantum numbers. In the following we characterize states 
solely by these two quantum numbers: the total component, S^, and the parity, P. An 
additional symmetry, corresponding to a global spin flip interchanging -1 and 1 while leaving 
alone, introduces another parity in the 5*^ = sector. Since most of our calculations have 
been performed in the = 1 sector we have not considered this symmetry. 

In Section |n| we shall briefly review the free boson model. Section |T| is concerned with 
our results for the energy gaps. In Section |I^ we discuss the bulk correlation functions. 
Section |V] addresses the effects that the open ends have on the system and discusses the 
boundary correlation functions. Theoretical and numerical results are presented for the two 
types of correlation functions in these two sections. In Section ^ we discuss the equal 
time structure factor in conjunction with experimental results and theoretical estimates. 
Section [VIl] discusses the single mode approximation for this problem. 



II. FREE BOSON MODEL 



It is believed that the isotropic S = 1 Heisenberg antiferromagnetic chain can be approx- 
imately described, at low energies, by the non-linear a (NLcr) mo delS, with Hamiltonian 



H 



^ / dx 



9 \ox , 
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2 J as, 



(2.1) 



where a is the lattice spacing. Since we have set ^ = 1, the velocity, v, has dimension 
energy times length. The massive triplet of fields (f) is restricted to have unit magnitude, 
(pi^ = 1, and and 1 describe the sublattice and uniform magnetization respectively. This 



Hamiltonian can be arrived at directly from a lattice HamiltonianEj, by use of the relations 
0(2i+l/2) = [S2i+i-S2i]/(2s), l(2i+l/2) = [S2i + S2i+i]/(2a). We have then approximately 
Si = (— 1)*S(/) + al, with 1=1/ {vg)4> x [d(f)/dt). The NLcr model has Lagrangian density 



1 

2^ 



1 fd^' 



— V 




(2.2) 



It is possible to obtain a linear model in the spirit of a phenomenological Landau- Ginzburg 
model describing the correct physics at the mean field level by introducing explicitly a mass 
term and a quartic term for stabilityfey while lifting the non-linear constriction (fp' = 1. The 
D term in Eq. (|1.1|) describing the single ion anisotropy will split the massive triplet into 
a low-lying doublet, with gap A^, and a higher lying singlet, with gap A||. One usually 
takes the z axis to be along the chain, describing the out of plane 1 1 mode, with the x, y 
axis perpendicular to the chain describing the in-plane _L mode. The phenomenological 
Lagrangian then becomes 



i=l 





2v- 



(2.3) 
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Here we have Sj ~ Sy^(— + al, with g = g±,g\\ for the perpendicular and parallel 
modes, respectively, and 

1 = i0 X ^. (2.4) 

V at 

where (p has been rescaled by a factor of Here we shall want to allow for the velocities 
for the different modes to be different in which case we have to generalize the above equation 
for 1 to obtain 

1- 1- 1- 1- 1- 1- 

1 = { — <Py4>z 02/02, — 020X- 020X, —(px4>y 4>x4>y)- (2.5) 

Vz Vy Vy 

A simple mean field theory in the spirit of the Landau- Ginzburg model can then be obtained 
by considering the free model with a = 0. This we shall refer to as the free boson model. 

The DMRG method that we shall use applies best to systems with open boundaries. 
The open boundaries has the effect of leaving a = 1/2 degree of freedom at each end of 
the chain, see Fig. |1], which will interact with the rest of the system. The open boundaries 
will thus introduce an interaction of the following form 

Hi = -X[s^(l>{l) ■ S; + al(l) ■ S; - s^(t^{L) ■ S'^ + al(L) ■ S'J, (2.6) 

where S'l and S'l^ are two S = 1/2 excitations known to exist at the end of the open 
chain0. Here we have implicitly assumed that the coupling to the staggered and uniform 
magnetization can be described by one coupling A, we shall, however, take A to be different 
for the in-plane and out of plane couplings, i.e. A^ 7^ A|| and equivalently g± 7^ g\\. In the 
following we shall therefore explicitly write ga, Xa where necessary. The individual spins, Sj, 
can now be represented as 

Si ^ s^i-iy-'cp + al + <5i,iS; + 6,^lS'l, (2.7) 

with 1 as in Eq. ( p.5|) . In the following we sometimes set the lattice spacing, a = 1. An 
alternative field theory treatment of the S = 1 Heisenberg chain using Majorana fermions 
has also been proposed^. This model predicts rather complicated forms for the correlation 
functions and we shall therefore use the conceptually simpler free boson model. 



III. ENERGY GAPS 

Due to the presence of the D term in Eq. ( |1 . 1| ) the Haldane gap will be split into a low- 
lying doublet, A_|_, and a higher-lying singlet, A||. These two magnon modes have quantum 
numbers of 

|a>=Ee^''^x|0>' (3-1) 

X 

where a = z for the singlet, and a = x or y for the doublet and |0 > is the singlet ground- 
state. Thus the singlet, \z >, has total 5"^ component, = 0, whereas the doublet, 
|± >= {\x > ±i\y >)/-\/2, has = ±1. The chains that we consider here are subject 
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to open boundary conditions and we shall choose the exponentially low-lying 1~ state as 
our reference state. In order to determine the two gaps we shall use a procedure similar to 
what was recently used for the isotropic chain@. The addition of a bulk magnon changes 
the parityS. This can be seen as follows. We only consider chains of even length. For 
these chains the ground-state is a singlet with even parity, 0^. Above the ground-state is an 
exponentially low-lying triplet, 1~. In the thermodynamic limit the triplet and the singlet 
become degenerate and the ground-state four-fold degenerate. This spectrum can be seen to 
arise from the two S = 1/2 end-excitations forming either an odd parity singlet or an even 
parity triplet, in addition to an overall parity flip coming from the rest of the ground-state. 
This parity-flip can be understood from the valence bond solid state! where we draw two 
valence bonds emanating from each site. These valence bonds represent singlet contractions 
of pairs of 5* = 1/2's so they have a directionality associated with them. When we make a 
parity transformation we flip the orientation of an odd number of valence bonds resulting 
in a (— ) sign. This is schematically depicted in Fig. |^. Thus, the parity, Pe, of a state with 
no magnons present is (+) if the end-excitations combine into the singlet and (— ) for the 
triplet. The parity of higher excited states, containing one or more magnons, is a product 
of three factors, PEPsv/Pm- Pm contains a contribution of (— ) from each magnon present. 
This is because the magnons are created and annihilated by the staggered magnetization 
operator, and this changes sign upon switching even and odd sublattices. Psw is the parity 
of the spatial wave-functions of the magnons. The gap to the doublet should therefore be 
calculated with respect to a state that has Sf. = ±1 and parity (—1) as compared to the 
reference state 1~. One candidate is thus the state 2"^. As shown in Ref. we expect the 
L dependence of this gap to be 

Our results are shown in Fig. |^ where the solid line indicates the best fit to the above form 
0.2998(1) + 105.2(1) (L - 1)"^ _ 794(2) (L - 1)"^ From this we can estimate A± and 

Ax = 0.2998(1), = 2.53(1). (3.3) 

The value of this gap previously obtained by Golinelli ahil of = 0.301 is in good agree- 
ment with the above result. v± can be compared to the value of v± ~ 2.53 that can be 
extracted from exact diagonalization results! of the energy as a function of k with D = 0.2. 
Presumably the dependence of v on D/ J is fairly small. 

The gap to the heavy magnon, A||, is more difficult to obtain. Following the above 
reasoning we use one of the excited states in the l"*" sector to calculate A||. The expected 
dependence on L is 

We show our results in Fig. |^ where the solid line indicates the best fit to the above form 
0.6565(5) + 42.5(l)(Li)-2 + 521(2)(L - l)-^. We can now estimate A|| and 

A|| = 0.6565(5), v\\ = 2.38(1). (3.5) 
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The value of A|| is in good agreement with previous resultsH^I. Extracting v\\ from exact 
diagonahzation resultsS with D/J = 0.2 we obtain f|| ~ 2.40 which supports our finding of 
W|| = 2.38(1) for D/J = 0.18. 

For periodic systems the energies will approach their asymptotic value exponentially fast. 
It is then easy to boost results for small chain lengths, to obtain rather good estimates of 
the asymptotic value, bv eliminating the exponential corrections through the application of 
Shanks transformationsE3. For the open ended chains, that we consider here, the corrections 
are powers in L to leading order given by Eq. (|3.2|) . In this case it is also possible to 
apply somewhat more complicated algorithms to extract asymptotic values for the gaps. 
We applied the alternating e-algorithm0 to our data and obtained values for the gaps in 
complete agreement with the above quoted values obtained by fitting to Eq. ( p.2|) . In 
general these techniques seems to yield a somewhat smaller accuracy for open boundary 
conditions&ii as compared to periodic systems. 

Using the DMRG method it is quite easy to extract the ground-state energy per spin 
in the thermodynamic limit by considering the quantityiiS (El - ^l-4)/4. We obtain 
eo/J= 1.2856861. 



IV. BULK CORRELATIONS 

We now turn to a discussion of the bulk correlations, < S^Sj >, where both i and j are 
far away from the boundaries. For the S = 1 systems we find using Eq. (|2.7|) 



< sfs^ >= i-iy+^ga < (j)''ixi)(i)''{xj) > W < r{xi)r{xj) > . (4.1) 



Here < ■ > denotes ground-state expectation values. Let us start with the first term. We 
shall follow the approach of Ref. We expand the staggered magnetization field, 4>, using 

4 

these operators 



a mode-expansion in the magnon operators afc,a[.. We use the relativistic normalization of 



al:^] = 47TVaUJkS{k-k'), (4.2) 



where Uk = + {vakf. With the definitions X = (Xo,Xi) = {t,x/v), K = (Ko,Ki) = 
(u;, vk) and K ■ X = cut — A;x, we write for the mode-expansion 

'A(x,t)=/^(e-^-^a, + e^^-M)- (4-3) 



Using the mode expansion together with the commutation relations Eq. ( [4. 2] ) we obtain 

r dk 1 

ga < O|0'^(a;, 0)0^(0,0)10 >= QaVa / — e-'^- , (4.4) 

J 47r y^A^T^ 

The integral can be expressed in a simple form using the modified Bessel function Kq, and 
we obtain 

Qa < O|0"(x, 0)0^(0,0)10 >= {x/ia) I^Po. ^"'^ (4.5) 

2^ 2J2T:\x\/ia 
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where we have set = Va/Aa. 

Due to the factor of 2 in the size of the correlation lengths only the leading Bessel 
function behavior, given by Kq, can be observed for the perpendicular bulk correlation 
functions, < S^Sj > and < SfSj >, as we shall see below. For < S^Sj > we find by fitting 
to the expected form 

< S^S^ >= (-l)^0.1843(2)iro (x/8.345(8)) , (4.6) 

with = 2.81. The results are shown in Fig. ^ where the solid line is Eq. ( [4.6|) . This fit 
determines the correlation length C,±, 

^± = 8.345(8). (4.7) 

The numerical values for this correlation function is listed in the appendix. In Fig. ^ it is 
also clear that the above form fails when the chain end is approached at the right hand 
side of the figure. We see that the relation = v±/A± is obeyed to within 1% if one uses 
the values for v± and Aj^ obtained in section |IT1[ An estimate of g± can also be extracted, 
gj_ ~ 1.16. Note that the Lorentz invariant form of Eq. ( [4.6|) works even for x < C,± where 
the asymptotic form of Eq. ( ^4.5| ) is not valid. 

The remaining term, < l°'{xi)l"'{xj) > constitutes minute corrections to the overall 
form given by the Bessel function Kq. P{x, t) contains four terms with two-magnon creation 
and annihilation operators. The only contribution to < l^{xi)l^{xj) > in the ground-state 
will come from the double creation term, l^. Thus the term < l°'{xi)l°'{xj) > in the 
correlation function effectively corresponds to two magnon excitations. From the mode 
expansion we see that 

i:ix, t)=^[ f'^^" (u,, - u;,.)e^^^'^^">''alUll (4.8) 

In this equation v± = = Vy. We then see, through the use of the commutation relations, 
that 

r rlh'rih" 

a2<0|r(x,0)r(0,0)|0>=a2 / -^p^(c.,,-c.,,OV(^'+^'>. (4.9) 



Rearranging the integral we obtain 



a" 



(4.10) 



With our definitions the last of the two integrals is equal to Ko{xA±/v±)/{7!'v±). The first 
of the two integrals can be evaluated through the relations 

r riy A 2 r A 2 ^2 

/ '^u,,e^''^ = ^ duVTT^^e-/^- = ^(-f^ + l)i^o(x/a), (4.11) 
J Ztc v±27r J v±TT ax^ 

where we have used the relation ^_|_ = v±/A±. By applying the well known recurrence 
relations for the modified Bessel functions it is easy to show that 
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l)Mx/U) = --l—K,{x/i^), (4.12) 



where Ki{x) = —dKQ/dx. Putting our results together we find that 



< 0|r(x,0)r(0,0)|0 >= -^r^iro(a;/a)^^^i(a;/a), x^O. (4.13) 

One should note that Ki is positive definite and the contribution is thus negative. This term 
will asymptotically behave as ~ — (1/x^) exp(— 2a;/^j_). Since 2/^± ~ ^|| the exponential 
dependence is of the same order as the leading term given by Kq. We have evaluated this 
term in the free boson model neglecting interaction effects, these could be of importance but 
are unlikely to change the asymptotic behavior of the two magnon term (see also section |V|). 
For the < SfSJ > correlation function this last term has a more complicated form with an 
exponential decay faster than the leading term and it is therefore more difficult to observe. 
Using these results we now consider the < S^Sj > correlations. From the above we have 

< SfS] >^ (-l)'^^o (^/^ll) - (x/U) ix/U) . (4.14) 

Here x = \i — j\, A = g\\, and B = 1 from the above derivation. Since is roughly twice ^|| 
in NENP there is a chance the last term is detectable in the < SfS^ > correlation function. 
In Fig. g is shown < S^qS? > as a function of |50 — ^l. Fitting to the above form we find 



< StS'j >= (-1)^0.2178(6)7^0 (x/3.69(l)) 



-Om7(3)Ko (x/8.9(9)) -^K^ (x/8.9(9)) , (4.15) 

\x\ 



with = 6.25. The solid line in Fig. ^ connects the discrete points obtained from Eq. ( [4. 151) 



The numerical values for this correlation function is listed in the appendix. This fit deter- 
mines , 

en = 3.69(1). (4.16) 

The relation = f||/A|| is then satisfied to within 2% if one uses the values for f|| and A|| 
obtained in section |ITT[ In addition we also find g\\ ~ 1.37. We can also compare the fitted 
coefficient for the two-magnon term, 0.007, with the theoretical prediction of 0^/(27?^^^). 
Setting a = 1 and using the previously obtained estimate for we see that the prediction is 
~ 0.006 in good agreement with the fitted value. The value we find for ^± is also in reasonable 
agreement with what we previously obtained by looking directly at the < SfSj > correlation 
function. 

In the case of the isotropic chain, where D = 0.0, we also expect the bulk correlation 



functions to behave as Eq. ( ^4.14| ). Our results are shown in Fig. ^ for < Sj^S^ >. Since 



the calculation is performed in the 1~ state we expect a disconnected term of the form 
< 5*70 >< 5"/ > which we subtract from < S^^S^ >. Such a disconnected term is also 
present when D = 0.18 but is orders of magnitude smaller due to the smaller correlation 
length and due to the fact that 5*50 is right at the middle of the chain whereas S^q is much 
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closer to the boundary for a 100 site chain. For small arguments in < S'|S'f > we expect the 
same effect to be present in the O"*" state since the wave function locally looks the same; only 
by considering large arguments will one be able to differ between the two states. Fitting 
< S^qS^ > - < S^qX > to Eq. §11) we obtain 

< S^S] > - <St >< >= (-1)^0.1999(1)^0 (a;/6.03(l)) 

-0.008(5)iro {x/7{2)) -^K, (x/7(2)) , (4.17) 

with a = 20. The fit works very well even for rather small values of \i — 70|. From this 
fit we can determine ^ = 6.03(1) (in complete agreement with Ref. and g = 1.26. The 
numerical values for this correlation function is listed in the appendix. For the isotropic 
Heisenberg chain the best numerical estimates for A, v, ^ ar A = 0.41050(2), V = 
2.49(1), ^ = 6.03(1), which shows that the relation ^ ~ f/A is correct to within 0.5%. 
Again we can compare the coefficient in front of the double Bessel function term to what 
we expect from the free boson theory, l/(27r^^) ~ 0.008, in very good agreement with the 
fitted value. The double Bessel function term is in this case rather difficult to fit due to 
the fact that it is exponentially small compared to the leading Kq term. The value we 
find for the correlation length from this term has therefore rather large error-bars as do the 
coefficient in front of it. Both are, however, in good agreement with the free boson theory. 
The leading Bessel function behavior, -KoiOf the bulk correlation functions has previously 
been determined for the isotropic casei'lrLl. 

The values g± ~ 1.16 and g\\ ~ 1.37 obtained above gives an average g = 1.27 which 
compares well with what we obtain for the isotropic chain, g = 1.26. The value of g for 
the isotropic chain can roughly be compared to the value of (7 ~ 1.44 for the bare coupling 
obtained!^ for the isotropic chain by 1/5* expansion. 



V. END EFFECTS 



As mentioned above, all of our numerical results are for chains subject to open boundary 
conditions. We therefore need to consider chain end effects. While at first sight such effects 
might seem undesirable they allow consistency checks on our results. For instance, they allow 
another way of estimating ^_|_ by looking directly at the energy levels. As mentioned, the 1~ 
state is exponentially low-lying. The gap between the 1~ state and the true ground-state, 
0"*", decays exponentially with the chain length, L. The energy difference arises because the 
two S = 1/2 end-excitations are in the triplet configuration in the 1~ state as opposed to a 
singlet for the ground-state, O"*". The interaction between the S* = 1/2 end-excitations are 
predominantly determined by the exchange of virtual magnons of the Aj^ typeEl. The energy 
difference should therefore decay exponentially with a characteristic length equal to C,±- This 
can be seen by integrating out (f> in the quadratic Hamiltonian of Eq. ( |2.3| ), including the 
interaction term between the chain-end spins and the staggered magnetization field, but 
ignoring their coupling to I for simplicity. (This leads to effects which decay more rapidly 
with L.) We obtain an effective action 5*05(8'-^, Sl)- Ignoring the interactions involving 1 in 
Eq. ( |2.6D we can write 
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We then find 

A 



5eff = -^9 J dndr,[S[{n) - S',{n)] . [S[{r,) - S',{r,)] 

dkdn V 



(27r)2 /t2 + v2f.2 + ^2 



where k is the imaginary frequency. By Taylor expanding S'^ l{t2) around T2 ~ Ti, we obtain 
a series of instantaneous terms involving increasing numbers of r derivatives. The lowest 
order term gives, after doing one r integral, an effective Hamiltonian: 

//.^A^A-Sl/gj^e'^'-". (5.3) 
The resulting integral now yields the familiar exponential form and we can write, 

H^s = \'gS[ ■ S'z.^^, (5.4) 

where ^, A and g are different for the two modes. Note that this effective Hamiltonian is 
exponentially small, for large L. It follows that a time derivative of S'^^ is of order this expo- 
nentially small energy scale. Hence the higher derivative terms in the effective Hamiltonian 
will be suppressed by powers of the ratio of this exponentially small energy to A and can be 
safely ignored. Due to the large difference between ^_|_ and ^||, we can write approximately 



1 r n P^iL-l)/^± 



Ai^?±^ S[^S', + S[-S'^\ . (5.5) 



2 

The equivalent contribution from the heavy magnon is smaller by a factor exp(— (L — 
l)(l/^ll — l/^±)), and can therefore be ignored. We now want to compare the 0"*" state 
with the 1~ state. In the 0+ state the wave function for the two spin-1/2 end excitations is 
approximately (l/v^)(| ti> — | iT>), while in the 1~ state it is just | tT>- We thus obtain 

In Fig. ^ we show this energy gap on a log scale as a function of the chain length L. Clearly 
there is an exponential dependence and we estimate by a simple linear fit to the data 
obtaining ~ 8.38(4) in nice agreement with the above result. The full functional form is 
given by Ei- — Eq+ ^ 0.27exp(— (L — l)/8.38). Using the above obtained values for Aj_ 
and g± we find A_l ~ 0.53. 

An equivalent analysis can be done for the isotropic chain. In that case we obtain 

3A2 e-(^-i)/« 

E,^-Eo. = —g^^. (5.7) 
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Fitting to the numerically determined energies we find Ei~ —Eq+ ~ 0.592 exp(— (L— 1)/6.07). 
The observed correlation length is again very close to what we found above using the bulk 
correlation function. The prefactor allows us to determine A. Using the values g ~ 1.26 and 



A = 0.4107, determined above, we find A ~ 0.72. 

We now consider < Sj > in the 1~ state. < S" > must vanish ioi a = x or y by rotational 
symmetry, but can be non-zero for a = z, having the value: 

< St >= 6n < > +6,L < S'l > +{-iy-'^ < >+a< F{x,) > . (5.8) 

We expect S'l ~ S'£ ~ 1/2 in the 1" state, from Hes of Eq. (|5l^). Due to the very small 
energy scale in ifefr the effective chain-end spins are "frozen" on time scales of 0(A~^) so we 
can treat them classically in considering the response of the bulk fields, 4> and 1. They act 
as classical sources at the two ends of the chain. If we temporarily ignore the couplings of 
the chain-end excitations to the uniform magnetization density, 1, then we can calculate the 
response of exactly, from the free boson model of Eq. ( |2.3| ) and ( p.6| ). The field equation 
for (h^ixj) is 



[^0^2 -^)^' = -^[^(^ - 1) - - ^)]- (5-9) 
This equation has the solution 



(2A|| 



(5.10) 



and we see that \ < Sf > \ should decay exponentially away from the chain ends with a 
correlation length ^|| = f||/A||. Note that, unlike the bulk correlation function, there is no 

1 / -y/jxj prefactor in this case; rather than a Bessel function, we obtain a pure exponential. 

Next we consider the uniform part, a < (xi) > . It is presumably not possible to find a 
closed form expression for this quantity, from the free boson model of Eq. ( ^.3| ) and ( |2.6| ) even 
after replacing the chain-end spins by their expectation values. We will content ourselves 
with doing lowest order perturbation theory in the coupling. A, between and S'^ to make 
a rough estimate of its magnitude. This gives: 

a < r{xi) >=^ jdT< P(x„ 0)P(0, r) >o, (5.11) 

where this expectation value is calculated in the theory with A = 0. This is similar to 
what was calculated in Sec. |^ on bulk correlations in Eq. ([4.13|) . Using the fact that 



Ko{yJx'^ + i;2r70 obeys the relation {-d'^/dx^ - d'^/dr'^ + l)KQ{y/x^Tv^ / = 6{x)6{t) 
we use Eq. (^l3| ) to obtain 

a < r{x,) >= £^ J dTKo{^^+^y^^)^Ko{^^^+^yu). (5.12) 

For large x the Bessel function has the behavior {1/ ^\x\/^±) exp(— — v^t"^ /{2\x\^±)), 
which results in an asymptotic decay for a < l^{xi) > as —{1/ x^^'^)e~'^^^^^^^ . 
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In Fig. P we show our results for | < 5*/ > | as a function of |1 — i\. The sohd hne 
connects the discrete points obtained by a fit to the form 

2{i-l)/8.45{6) 

(-l)'-i0.380(2)e-(*-i)/3-^°3(^) + 0.133(4)^ ^= , (5.13) 

with = 0.06. Again we see that the value for the correlation lengths, 

^11 = 3.703 ± 0.004, ^± = 8.45, (5.14) 

are in excellent agreement with what we determined from the bulk correlation functions. 
The absence of any 1/x prefactor in the staggered part is also in accord with the boson 
model. However, the l/\/x prefactor in the uniform part is in disagreement with the 
prediction of the boson model, also the sign is wrong for that term. We do not understand 
these discrepancies. We cannot exclude the possibility that they are somehow generated by 
the iterative numerical procedure, although we find it unlikely. From the staggered part of 
the fit we can determine A|| ~ 0.85. 

In Fig: ^ we show < Sf > for a 100 site isotropic chain. In this case we fit to an 
exponential form and we obtain with a = 6.7, 

0.486(1) [exp(-(i - l)/6.028(3)) - exp(-(100 - i)/6.028(3)]. (5.15) 

One should here remember that the two contributions from each end of the chain add up out 
of phase, hence the minus sign between the two exponential terms. This clearly establishes 
^ = 6.03(1), again in complete agreement with our previous result and Ref. |^. Also we 
can estimate A to be A ~ 0.71, which is in excellent agreement with the value A ~ 0.72 
determined above. 

Now we consider the correlation function, < S^S^ >, with one operator at the edge of the 
system. This has a large number of contributions using the expression of Eq. (|2.7|) for the 
spin operators and the Hamiltonian of Eq. (|2.3| ) and (|2.6| ). Let us first consider < SfSf >, 
which is a bit simpler because < >= 0. The staggered part is given by: 

i~iy~V < 0"(xi)0"(x,) > < s[^rixi) >]. (5.16) 

While the first term gives a similar expression to that obtained in the bulk, the second term 
has a somewhat different asymptotic behavior. This can be extracted by doing first order 
perturbation theory in A: 

< S[^r{xi) >~ Kg. I dr < ril,r)r{x,,0) >o< 5f (r)5r(0) >o • (5.17) 

Here the Green's functions are calculated in the A = limit. Hence, S'^(r) is time- 
independent and < 5f (r)5f(0) >o= 1/4. Thus: 

< J 2^ .iP + Ai = — 2A^- 

This has a pure exponential decay, unlike the bulk part which has a 1/y/x prefactor. The 
uniform part has an exponential decay of the form e~^^^^^^-^^^^^^^\ This decays much more 
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rapidly than the staggered part and is essentially unobservable in our numerical work. Our 
data for < SfSj > are shown in Fig. 10. For small arguments there is an indication of a 



uniform part but we are not able to determine this. For large arguments the points of this 
boundary correlation function clearly falls onto one line and we make a fit of the following 
form 

< S'fS^ >~ (-l)^[0.277(2)e-^'/s-3^(2) ^ o.064(3)iro(a;/8.37(2))], (5.19) 

with X = \i — 1\, which is shown as the solid line in Fig. |1^. For this fit = 32. This fit is 
just marginally acceptable hinting at the presence of a uniform part. We can again read off 
the correlation length 

a = 8.37(3), (5.20) 

in excellent agreement with what we previously obtained. The coefficient in front of the 
exponential should be ~ 0.26 using our previously obtained values which is close to 
the fitted value. 

The boundary correlation functions < S^S^ > is more complicated for two reasons. 
Firstly, < S- >^ a.s discussed above. Secondly, the uniform part decays with the exponent 
2x/^± which is actually slightly sma//er than the one occurring in the staggered part, x/^\\. 
Therefore it is also observable. We consider the connected part: < S^S^ > — < Si >< >. 
The staggered part is expected to have a pure exponential term as well as a Bessel function 
term by similar arguments to those given above for the boundary correlation function < 
SfS^ >. The uniform part is expected to have a l/x^^^ prefactor, by arguments similar to 
those given for < Sf >. In Fig. |TT| we show our results for < SfS^ > — < >< >. As 
is clearly seen in this figure this correlation function not only has a staggered part but also 
a uniform part. We obtain 

2(i-l)/8.33(3) 

< S^S^ > - <Sl >< >~ (-l)'-^0.117(l)iro((^ - 1)/3.66(1)) - 0.081(1) ^== — , 

Vi — 1 

(5.21) 

with = 0.31. Our data are clearly consistent with both the Kq term and a uniform term 
being present in < S^S^ >. Furthermore the value we obtain for as well as for ,^||, 

eil = 3.66(1), (5.22) 

is in very good agreement with what we obtain for the same correlation length from the 
bulk correlation function, < S-Sj >. However, there are unexplained discrepancies with 
the predictions of the free boson model. A pure exponential alternating piece was not 
measured and the uniform part decays with a I/a/x prefactor, rather than 1/x^^'^. We 
cannot exclude the possibility that this discrepancy somehow is generated by the iterative 
numerical procedure, although we find it unlikely. 



In Fig |T2| is shown the boundary correlation function < S^S^ > — < >< > for a 



100 site isotropic chain. In this case we fit to the simple form 

< S^S^ > - <Sl >< St >~ (-l)^-^0.07(l)iro((z - l)/5.93(3)), (5.23) 
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with = 13.7. In this case we do not see an exponential term in the staggered part nor 
do we find any uniform part. The obtain correlation length is in reasonable agreement with 
what we obtained from the bulk correlations. In this case we were not able to fit the uniform 
part of the correlation function. 



VI. EQUAL TIME STRUCTURE FACTOR 

Although, as mentioned, it is not possible to obtain direct information about dynamical 
properties using the DMRG method, it is in principle possible to obtain the equal time 
structure factor, by simply Fourier transforming the bulk correlation functions. Exact diag- 
onalization results for S{k) for shorter chains, L = 16, allows us to check the validity of this 
procedure. Here we use the standard definitioni 

S{k, u;) = — V / dte^'^'e-''''' < S{r, t) ■ S{0, 0) > . (6.1) 
2% ^ J 

With this definition the equal time structure factor can be written in the following way 

S{k) = j dujS{k, uj) = J2 < 5(r, 0) ■ ^(0, 0) > . (6.2) 

In particular, we use for our data the relation 

S{k) =< S{50, 0)S(50, 0) > +2 ^ g-ife(50-r) ^ g^^^ 0)5(50, 0) > . (6.3) 

r 

In Fig. |T3|, |14| this is done so as to obtain 5*" (A;) and S^{k), respectively. Our results 
are shown as the open squares and circles, respectively. The real space bulk correlation 
functions. Fig. ^, ^ we have used in obtaining the structure factors were calculated for an 
open ended chain and one might wonder about the validity of this approach. However, since 
we are able to obtain results for very long chains the approximation made by assuming that 
the correlation functions are equal to the correlation functions for a periodic chain with the 
same length is very small. For 5*11 the error is negligible since the associated correlation 
function decays very fast. However, for S-^, as can be seen in Fig. there are some end 
effects present in the correlation function which we estimate to give rise to an error of 0.05% 
in S^{k = vr), by using the fitted Bessel function form instead of the actual data points 
for r < 20. Away from k = ir we expect the error to be of the same order of magnitude, 
although somewhat larger near k = 0. The advantage of using the simple Fourier transform 
is that the full weight of the structure factor is conserved. That is the total moment sum 
rule, 

J E[25^(A:) + SHk)] = SiS + 1) = 2, (6.4) 
^ k 

is obeyed to within numerical precision. 



The dotted lines in Fig. |T3|, |T4| are square root Lorentzians (SRL). We define this as 
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9a^ I ^ , , , . ia = Va//^at (6.5) 

2 ,^l + (A;-7r)2(e,)2 

which is just the Fourier transform of the Bessel function form of the staggered part of the 
correlation function, Eq. (|4.4| ). In both figures we have used the values previously obtained 
for Qa and ^a- From section |I^ we know that the SRL form should work very well close to 
k = n, since Qa, C,a was obtained by fitting to the staggered part of the correlation functions. 
From Fig. |13|, |l^ it is, however, apparent that the simple SRL form only describes the 
structure factor well from fc/7r~0.7tofc/7r = l. 

For the isotropic chain the ground-state is a singlet and therefore J2i S{i,t)\0 >= 0. 
Thus the structure factor must vanish as k"^ close to A; = 0. If we introduce anisotropy then 
S*" still vanishes as fc^ close to A; = but S-^ can take on a non zero value at A; = 00. 
Within the frame work of the free boson approximation it is possible to derive approximate 
expressions for the structure factors in the presence of anisotropy. This has previously been 
done assuming uniform velocities v = Vx = Vy = Vz i^n which case simple expressions for 
S{k, uj) can be derivedS. We now generalize these results to allow for different velocities but 
focus only on the equal time structure factor, S{k), since expressions for the full dynamical 
structure factor, S{k, u), in this case will be rather complicated. The structure factor close to 
k = n will be completely dominated by correlations in the uniform part of the magnetization 
density. If we use the general expression for I, Eq. (\i.5i\), along with the mode expansion, 
Eq. ([4.3|) , and the commutation relations, Eq. ( [4.2|) , we find that I will contain a double 
creation term of the following form for (compare Eq. (|4.8|)) 



i:{xtt) = . / J^Mfc. ^^iK^^K.yx^^ _ 5^)4 at (6.6) 

Again this term will be the only term contributing to < l^{xi)l^{xj) > in the ground state. 
We thus obtain for the equal time correlation function, 

< o|r(x,o)r(o,o)|o >= / -JM!^e'iky+k.)^(^^ _ ^)\v^. (6.7) 

J IQ-K^UJky^k, Vz Vy 

Thus we see that the two magnon contribution to S-^{k) = S^^{k) = S^^lk) is 
2 f dkydkz ,Uk,Vy uJkyVz 



J StT UkyVz UJk.Vy 

For the parallel mode the two magnon contribution to the structure factor, ^"(/c) = S^^{k), 
becomes simply 

2 [ dk^dky UJk^ ^ky 



a' / :^^{^ + ^ - 2)5ik -kx- ky), (6.9) 



in LUky iOk 



which is just the Fourier transform of Eq. (|4.13| ). Here the overall scale is not a free param- 
eter, as opposed to Eq. ( |6.5| ) where ga is a free parameter, because / dxl = must obey 
spin commutation relations. These two predictions are shown as the solid lines in Fig. 
Qualitatively the free boson estimates seems to agree well with the numerical results for k 
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close to 0. For fc/vr > 0.1 the free boson estimate for the two magnon part of S{k) is larger 
than the numerical results for the full structure factor. Since the free boson theory doesn't 
take interaction effects into account this is perhaps not too surprising. 

In Fig. |13|, |14| are also shown the inelastic neutron scattering (INS) results of Ref. p!7 
which are directly comparable to our results. The open triangles are points where to within 
experimental accuracy S*" and S*-*" were identical. The full squares and circles are data 
points where S*" and S*-*", respectively, could be resolved experimentally. Good agreement 
is obtained with the experiment apart from an overall scale factor of about 1.25 which 
was to be expected since the experimental total intensity exceeded the exact sum rule: 
(1/-^) J2ak 'S'""(fc) = s(s + 1) by 30% (± 30%). A very nice agreement between the numerical 
results and experimental data is evident. 

Fig. |15] summarizes our results for S"" and S""*". In this figure the two structure factors are 
plotted together with the different theoretical estimates. For the range (0.1 -0.85)A;/7r 
is the larger of the two, and only when k/ir > 0.85 or fc/vr < 0.1 does S""*" become dominant. 
The crossing of the two structure factors near k/n ~ 0.85 is correctly described by the 
SRL if different velocities for the two modes are allowed for. The crossing near /c/tt ~ 0.1 
seems also to be predicted by the free boson theory estimate for the two magnon part of 
the structure factor. This seems also to be supported by exact diagonalization studies for 

L = lei. 

In the case of the isotropic chain, with D = 0.0, it is possible to obtain an exact expression 
for the two magnon part of S{kuj) within the frame work of the NLcx model0 thereby taking 
into account interaction effects. In Ref. ^ it is shown that with S = = S"^ = S^, the two 
magnon contribution to S is given by 



Here K ■ K = 4A^cosh^(^/2) and |G(6')p is given by the expression 



|^..^|. _ 1 + jOM' / tanh(g/2) y 

This calculation contains no free parameters. The expression for the two magnon contri- 
bution to S{k,uj) can now easily be integrated numerically over u to obtain S{k). Our 
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numerical results are shown in Fig. along with the SRL form using the previously ob- 
tained values for v, g and ^. The SRL form is shown as the short dashed line. Also shown, as 
the solid line, is the NLa model results for the two magnon part of S{k). As seen in Fig. ^ 
there is excellent agreement between the theoretical and numerical results. The long dashed 
line is the prediction from the free boson theory for the two magnon contribution to the 
structure factor. The inclusion of interaction effects changes the shape of S{k) and we see 
that the free boson estimate, although qualitatively correct near k = 0, somewhat overesti- 
mates S{k) for larger k. Presumably the NLa model also predicts a 4-magnon (and higher) 
contribution to S{k,uj). This is not known exactly and not included in Eq. (|6.10|) . Hence 
the full S{k) in the NLcx model should be somewhat larger than Eq. ( |6.10| ). The very precise 
agreement with the numerical results may indicate that this multi-magnon contribution is 
very small. Alternatively, it may indicate that the NLa model somewhat overestimates the 
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two-magnon part. Our numerical results for S{k) are also in good agreement with previous 
results using Monte Carlo techniqueSill for chains of length 64. 



VII. SINGLE MODE APPROXIMATION 

As already mentioned the DMRG method does not allow us to obtain information about 
the dynamical properties of the system. However, if S{k, u) is dominated by an intense single 
mode, as is the case for k close to tt, then we can obtain an approximate dispersion relation 
through the use of the single mode approximation (SMA). This has previously been done 
for the isotropic chain@. This can then be compared to exact diagonalization result sil'iifl 
and to QMC results for the isotropic chaiJii. We assume that for both modes 

S{k, u) = So{k)6{uj — Uk) + S{k, u), (7.1) 

where S{k,uj) is non zero only for u > Uc and Uc > ujk- Here uOc denotes the bottom of the 
continuous part of the spectrum. We then look at the integral 

— / duuS{k,uj) = So{k) H / dujujS{k,uj). (7.2) 

CUfc J —oo OJk J ~oo 

Since, by assumption, S{k,u) only is non zero when u; > tUc > ci^fc we must have 

— / dujLoS{k,Lo) > / dujS{k,uj). (7.3) 

We can then write 

1 poo roo _ 

— / dujiuSik,uj)>Soik)+ dujSik,uj) = Sik), (7.4) 

LUk J —oo J —oo 

where S{k), per definition, is the equal time structure factor. The first moment of the 
dynamical structure factor obeys a simple sum rule^'" 

c/cucj^""(A;, uj) = -l< [[H, 5"(A;)], S'^i-k)] >, (7.5) 

D 2 

with equivalent relations for S^^ik^uj) and S^'''{k,uj). Here we have used 

i 

Explicitly calculating the double commutator we get the following relation 

/oo 
duuoS''''{k, uo) = -J[Fy - cos{k) + - Fy cos{k) + D{G, - Gy)] 
-oo 

/oo 
duuS'^k, u) = -J[F^ - Fy cos(fc) + Fy - F^ cos(A;)], (7.7) 
-oo 

where Fa =< S^S^^-^ > and Ga =< (5'f )^ >. From Eq. ([7^ ) we then get the SMA dispersion 
relations 
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u^i < cojMAik) = -{J{Fy + - cos(fc)) + D{G. - Gy)}/S^ik) 

4 < 4MAik) = -{J{F, + Fy){l - cos{k))}/S\k). (7.8) 

We see that we always have ujswia > ^k, where Uk describes the dispersion of the singular 
mode. The inequality is only fulfilled as long as uj^ is below the edge of the continuous part 
of the spectrum, c^c, see Eq. ([T.!]). If we denote the lower edge of the spectrum by oj^, then it 
IS easy ^SMA > always is satisfied independent of the position 

of any singular modes. This simply follows from the trivial inequality 

— dujuS{k,uj)> dujS{k,u) = S{k). (7.9) 

Retracing the above steps we find ujsma > t^e- 

The long-distance field theory limit of the antiferromagnetic Heisenberg model is the 
0(3) NLo" model, Eq ( |2.2| ), describing a massive triplet of fields. The unit vector of fields 
describes the staggered magnetization of the spin chain. The spectrum for the S = 1 spin 
chain at /c = tt should therefore be dominated by large single modes corresponding to the 
massive triplet together with a small continuum from 3,5,. . . magnons. The NLa model has 
no bound states and the spectrum at A; = must therefore be a continuum corresponding to 
the excitation of two massive particles with k = ±7r together with a small contribution from 
4,6,. . .magnons. If anisotropy is present the spectrum at A; = tt consists of a light magnon 
with mass A_l and = ±1 and a heavy magnon with = and mass A||. At A; = the 
lower edge of the continuum of two magnon excitations for out of plane modes, cu", which 
have = with respect to the ground-state, must be given by the excitation of two light 
magnons with k = ±7r and = ±1, corresponding to a gap of 2A_i_. For in-plane modes, 
cj-*-, which have = ±1 with respect to the ground-state, the lower edge of the two magnon 
continuum at A; = must correspond to one light magnon with = ±1 and one heavy 
magnon = and opposite momentum. Thus the gap should in this case be A_|_ + A||. 
One expects that somewhere between A; = and k = ir the lowest energy excitations will 
change from the single mode excitations at A; = vr to the two magnon continuum close to 
A; = 0. This picture of the spectrum at A; = is in good agreement with exact diagonalization 
result Ji. 

For the anisotropic Heisenberg chain, with D/J = 0.18, we have = Fy = —0.4969462 
and Fz = —0.4035174. In conjunction with < (5*50)^ >~ 0.6206864 this agrees with the pre- 
viously obtained value for cq/ J = —1.2856861 to within numerical precision. Using Eq. ( [7.8[ ) 
along with < {S^qY >= 0.6896568 we can now obtain the dispersion relations within the 
SMA. Our results are shown as the dashed lines in Fig. along with the experimentally 
determined dispersion relation from Ref. 0. The long dashed line is and the short 

dashed line is Ugj^jj^. The INS data from Ref. |17| are shown as solid squares and circles for 
the data points where a;" and uj-^, respectively, could be resolved and as open triangles for 
the data points where to within experimental accuracy a;" and uj-^ were identical. The two 
previously determined velocity, vu and v±, gives us dispersion relations close to tt which are 



shown as solid lines. Also indicated in Fig. |17| is the lower edge of the spectrum at A; = 0, 



oj^ = 0.986J and ui]^ = 0.60J, and at A; = 7r/2, uij- = 2.75 and lj^J = 2.65J, obtained from 
the exact diagonalization of Golinelli at al0. The numerical data is scaled with J = 3.75 
meV. The agreement between the experimental results and our dispersion relations (solid 
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lines) close to /c = tt is excellent. The dispersion relations are also in very good agreement 
with the experimental results of Ref. |^. The solid lines indicate a crossing of a;" and oj-^ at 
approximately k/ir 0.8. This has also been observed in exact diagonalization studiesH 
and QMC studied^ and seems to agree with the SMA prediction as well as with the fact that 
experimentally uj-^ and tu" only can be resolved for kj-n > 0.9. As expected we see that losma 
consistently is larger than the experimental results and our relativistic dispersion relations 
when k is close to vr. Thus there is clearly a small but non zero contribution to Sik,^ lS) from 
multi magnon processes at /c ~ vr. '^\ma becomes smaller than the experimental points 
at around k/n 0.7. However, cu^^v/a coincides with the experimental data all the way to 
k/n = 0.3. One should note that a complete agreement between ojsma and the experimental 
data doesn't necessarily mean that the single mode approximation is exact but can be due 
to the fact that the singular mode is no longer below the continuous part of the spectrum. 
At fc/vr = 0.5 both SMA estimates are slightly above the lower edges u;^'" obtained from 
exact diagonalizationii. The experimental and numerical data seem consistent with a well 
resolved single mode all the way to fc/vr = 0.3, but this may no longer be well separated 
from a small multi magnon continuum. The fact that our results for ujsma remains rather 
close to the experimental data is consistent with a large singular contribution to S(k, uj) in 
this region. At A; = we can only compare to the exact diagonalization results of Ref. |2^. 
The SMA results are now at almost twice the value of these results consistent with a large 
contribution from multi magnon processes. It therefore seems likely that the continuous 
part of the spectrum will only set in at fc/vr < 0.3 in NENP. This is not inconsistent with 
our findings for the structure factor. 

For the isotropic chain we have = Fy = Fz = eo/3 = —0.4671613 again in very good 
agreement with the knownB value of eo/J = -1.401484038971(4). In Fig. |T8|the SMA results 
are shown as the dashed line. Also shown in Fig. |18|, as open squares, are the QMC results 
of Ref. 1^. The QMC results should determine the lower edge of the excitations spectrum 
regardless of whether this is a single mode or a continuum. The previously determined 
velocities and gapil gives a relativistic dispersion relation close to ~ vr which is shown 
as the solid line. The agreement between the relativistic dispersion relation and the QMC 
results close to A; ~ vr is excellent as one would expect. The SMA results indicate again 
a small but nonzero contribution to S{k, u) from the multi magnon continuum close to 
k = n. For a large region around k = n the SMA seems to work well as compared to 
the QMC results. However, close to = the SMA is considerably larger than the QMC 
results indicating a large contribution to S{k, lj) from multi magnon processes. Also shown 



in Fig. is the bottom of the two magnon continuum, SyA^ + v'^{k/2y, valid near k = 0. 
This result is in very good agreement with the QMC results. 



VIII. DISCUSSION 

Using the DMRG method it is difficult to get a good handle on what the actual error- 
bars are. Presumably the finite lattice method that most of our results have been obtained 
with effectively works as a variational methodiiS. The error-bars quoted in the previous 
sections are therefore only statistical error-bars obtained by estimating the errors occurring 
in the DMRG. This does not include any systematical errors that are present. It is also 
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not a trivial matter to fit data obtained on finite lattices to theoretical forms that are only 
asymptotically valid. We therefore here recapitulate our main findings with error-bars we 
believe are close to reality but more or less phenomenologically obtained by judging the 
variation of the various parameters between different fits. 

A||/J = 0.6565(5), t;||/J= 2.38(1), ^||/a = 3.69(5), ^|| ~ 1.37, A|| ~ 0.85 

Ai/J = 0.2998(1), t;x/J= 2.53(1), ^^/a = 8.35(7), ^ 1-16, ^ 0.53. (8.1) 

Here all the results were obtained for an anisotropic chain with D/J = 0.18. The bulk 
correlation functions show evidence of a small but non zero uniform part. Close to k = % 
the equal time structure factor is very well described by a square root lorentzian. The equal 
time structure factor S-^ does in this case no longer approach as /c — >■ 0. In this region 
both S"" and S-^ are well described by a two magnon contribution to the structure factor 
calculated within the free boson theory. The agreement between the free boson theory and 
the numerical results is however only qualitative and judging from the comparison between 
a single mode approximation and the experimental data it seems likely that a single mode 
persist all the way to fc/vr = 0.3. By comparing the single mode approximation to exact 
diagonalization results it is seen that between /c = and /c/tt = 0.3 the continuous part of 
the spectrum should become dominant. 

For the isotropic chain, with D/ J = 0.0, we findi 

A/J = 0.4107(1), t;/J = 2.49(1), ^||/a = 6.03(1), ^-1.26, A ~ 0.72. (8.2) 

Also in this case we find evidence for a small but non zero uniform part of the bulk correlation 
functions. Close to k = tt the equal time structure factor is very well described by a square 
root lorentzian. In the region close to A; = the equal time structure factor can be compared 
to an exact result obtained from the NLo" model for the two magnon contribution to the 
structure factor. An excellent agreement between theory and the numerical data is seen out 
to fairly large values of k. This seems consistent with the behavior of the spectrum. Here 
the single mode approximation as compared to quantum Monte Carlo results indicates the 
onset of a continuous part earlier than for the anisotropic chain. 

From the results presented here it seems likely that a contribution to the structure factor, 
near A; = 0, from two magnon excitations should be observable in NENP since S-^{k = 0) 
is non zero. Experimentally no indication of this has been observed so far. Regnault et 
alii have performed neutron scattering experiments near k = and find no indication of a 
nonvanishing structure factor aX k = 0. Ma et al0 analyze their INS data obtained in the 
range k/ix > 0.3 solely in terms of a single mode and do not observe the expected two magnon 
continuum. These authors point out that the free boson two-magnon prediction for S{k,u), 
when convoluted with the experimental resolution function, yields a much broader peak than 
what is observed experimentally, at the smallest value of k studied experimentally, k = .Sn, 
although the integrated intensity is roughly correct. In Fig. |l^ we compare the unconvoluted 
free boson two-magnon prediction at = .37r with the experimental result. Note that 
the integrated intensity and location of the intensity maximum are approximately correct. 
Also note that the width of the theoretical curve is about 2 meV whereas the experimental 
width, of about 1 meV, appears to be resolution limited. Under these circumstances, the 
possibility that the experimentally observed peak is of two-magnon type, with a width 
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somewhat reduced by interaction effects, cannot be ruled out. Including interaction effects 
in the NLcr model does lead to a narrowing of the two-magnon peak as shown in Ref. 

The much larger width of the free boson prediction, of about 5 meV, shown in Ref. ^7\, Fig. 
(Ic), is primarily a consequence of convolution with the experimental resolution function. 
Although this resolution is very narrow in frequency, it is rather broad in wave-vector, 
A A; ^ .2571. The relatively narrow convoluted width of about 1 meV, was obtained in 
the experiment by tilting the resolution ellipsoid so its axis was approximately parallel to 
the dispersion curve. Under these circumstances, the result of convoluting any theoretical 
prediction with the resolution function is extremely sensitive to the precise shape of the 
theoretical dispersion relation (ie. the peak center vs. k). It is clear from Figure 18 that the 
relativistic two-boson dispersion curve starts to fail badly for k > .3tt. Hence convoluting this 
function with a resolution function of width .257r produces an enormous width in frequency. 
This is especially so since the intensity increases rapidly with k. 

Therefore, the very broad resolution of the experiment prevents a definitive test of the 
theoretical model. The extremely good agreement with the NLa model two-magnon predic- 
tion for S{k) shown here, and the shifting of spectral weight above the lowest state, shown 



for L < 18 in Refs. p3|,|8| suggest that the excitations may start to exhibit a two-particle 
character at around k ^ .Svr. Further experimental work with higher resolution and lower k 
will be needed to settle this question. 
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TABLES 



TABLE I. 


The bulk correlation functions for a 100 site S = 


1 antiferromagnetic Heisenberg 


chain with single ion anisotropy Z)/J = 0.18. 
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TABLE II. The bulk correlation function for a 100 site isotropic S = 1 antiferromagnetic 
Heisenberg. 
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FIGURES 

FIG. 1. Schematic drawing of the valence bonds showing the two 5 = 1/2 end-excitations. 
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FIG. 2. The gap Aj^ as a function of 100/(L — 1)^. The soUd hne indicates a best fit to the 
theoretically expected form. 

FIG. 3. The gap A|| as a function of 100/(L — 1)^. The solid line indicates a best fit to the 
theoretically expected form. 

FIG. 4. The bulk correlation function, < SfS^ >, for a 100 site anisotropic chain. Shown is 

< S^^Sf > as a function of |50 — z|. The center of the chain is thus at the left hand side of the 
figure while the chain end is approached to the right. The solid line represents a fit to the Bessel 
function form described in the text. 

FIG. 5. The bulk correlation function, < SfS'j >, for a 100 site anisotropic chain. Shown is 

< S^QSf > as a function of |50 — i\. The center of the chain is thus at the left hand side of the 
figure while the chain end is approached to the right. The solid line connects the discrete points 
obtained by fitting to the form described in the text. 

FIG. 6. The bulk correlation function, < SfSj >, for a 100 site isotropic chain. Shown is 

< SjqS^ > — < SjQ >< Sf > as a function of |70 — i|. The chain end is approached to the right. 
The solid line connects the discrete points obtained by fitting to the form described in the text. 

FIG. 7. The gap, E^- — Eq+, as a function of chain length. The dependence is exponential 
defining a characteristic length of 8.38. 

FIG. 8. < S? >, for a 100 site anisotropic chain as a function of |1 — i\. The end of the chain 
is at the left hand side of the figure. The solid line connects the discrete points obtained from the 
fit (-1)^-10.380(2) exp(-(z - l)/3.703(4)) + 0.133(4) exp(-2(z - l)/8.45(6))/V«^- 

FIG. 9. < Sf >, for a 100 site isotropic chain as a function of |1 — i|. The end of the chain is 
at the left hand side of the figure. The solid line represents a fit to the points with of the form 
0.486(l)[exp(-(z - l)/6.028(3)) - exp(-(100 - i)/6.028(3)]. 

FIG. 10. The boundary correlation function, < SfSf >, for a 100 site anisotropic chain. Shown 
is < SfSf > as a function of |1 — The end of the chain is thus at the left hand side of the figure. 
The solid line represents a fit of the form 0.277(2) exp(-(i-l) /8.37(2))+0.064(3)ii:o ((z-l)/8.37(2)). 

FIG. 11. The boundary correlation function, < SfSf > - < Sf X Sf >, for a 100 site 
anisotropic chain. Shown is < SfSf > — < Sf >< Sf > as a function of |1 — i\- The end of the 
chain is thus at the left hand side of the figure. The solid line connects the discrete points obtained 
by fitting to the form described in the text. 
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FIG. 12. The boundary correlation function, < SfSf > - < >< >, for a 100 site 
isotropic chain. Shown is < SlSf > — < Sf >< Sf > as a function of |1 — i|. The end of the chain 
is thus at the left hand side of the figure. 



FIG. 13. The structure factor, as a function of fc/vr for a 100 site anisotropic chain. The 
numerical results are shown as open squares. INS data from Ref. 17 are shown as solid squares 
for the data points where experimentally could be resolved and as open triangles for the data 
points where to within experimental accuracy and were identical. The dashed line is the 
SRL form obtained by fitting close to = tt. The solid line is the prediction from the free boson 
theory for the contribution to S"" stemming from two magnon excitations. 



FIG. 14. The structure factor, S"-*-, as a function of fc/vr for a 100 site anisotropic chain. The 
numerical results are shown as open circles. INS data from Ref. 17 are shown as solid circles for the 
data points where experimentally could be resolved and as open triangles for the data points 
where to within experimental accuracy SW and were identical. The dashed line is the SRL form 
obtained by fitting close to A: = vr. The solid line is the prediction from the free boson theory for 
the contribution to S"" stemming from two magnon excitations. 



FIG. 15. The structure factors, and S-^ as a function of /c/vr for a 100 site anisotropic chain. 
The numerical results are shown as open squares and circles, respectively. The dashed lines are 
the SRL form obtained by fitting close to /c = vr. The solid lines are the predictions from the free 
boson theory for the contribution to the structure factor stemming from two magnon excitations. 



FIG. 16. The structure factor, S, as a function of /c/vr for a 100 site isotropic chain. The 
numerical results are shown as open squares. The short dashed line is the SRL form obtained 
by fitting close to A; = vr. The solid line is the exact prediction from the NLcr model for the 
contribution to the structure factor from two magnon excitations. The long dashed line is the free 
boson prediction for the two magnon contribution to the structure factor. 



FIG. 17. Dispersion of the excitations as a function of k/iT for a 100 site anisotropic 
chain. 



-{k) 



The points are INS data from Ref. 17 
jJl^l.+vlAk- 



The solid lines are the dispersion relations 
7r)2, for the two modes valid close to A: = vr. A value of J = 3.75meV 



has been used to scale the data with respect to the experimental results. The long dashed line is 
'^SMA' short dashed line is ijJsMA- "^^^ fo^'^ arrows indicate the lower edge of the excitation 

spectrum at /c = 0, uj^ = 0.986J and wi' = 0.60J, and at /c = 7r/2, u)-^ = 2.75J and u)e = 2.65J, 
from the exact diagonalization results of Ref. 23. 
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FIG. 18. Dispersion of excitations as a function of k/ir for a 100 site isotropic chain. The points 
are QMC data from Ref. ^. The sohd hne is the dispersion relation uj{k)/J = a/A'^ + v-^{k — tt)-^, 
valid close to /c = vr. The long dashed hne is the bottom of the two magnon continuum, 
2a/A'^ + v^{k/2)'^ , vahd near k = Q. The dashed hne is the SMA prediction for the dispersion 
relations obtained from the equal time structure factor. The poor agreement between the SMA 
prediction and the QMC results at small k indicates the two magnon nature of the excitations 
there. 

FIG. 19. S{k = O.Stt, w) as a function of to for the free boson model (solid line) and INS results 
(open triangles) from Ref. |l^. We have used J = 3.75meV. The dynamical structure factor is 
obtained as a sum 5 = (1 + cos e)S^ + sin^eS\\, with cos = 0.1789 in order to compare with the 
experimental data. 
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Fig 3 Sorensen, Affleck 
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Fig 8 Sorensen, Affleck 
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Fig 12 Sorensen, Affleck 




Fig 13 Sorensen, Affleck 




Fig 15 Sorensen, Affleck 




Fig 16 Sorensen, Affleck 
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Fig 17 Sorensen, Affleck 
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Fig 18 Sorensen, Affleck 
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Fig 19 Sorensen, Affleck 



